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Abstract. Species extinction occurs regularly and unavoidably in ecological systems. 

[-J*] The time scales for extinction can broadly vary and inform on the ecosystem's stability. 

n , We study the spatio-temporal extinction dynamics of a paradigmatic population model 

^ where three species exhibit cyclic competition. The cyclic dynamics reflects the non- 

•^H equilibrium nature of the species interactions. While previous work focusses on the 

1^ coarsening process as a mechanism that drives the system to extinction, we found that 

O^ unexpectedly the dynamics to extinction is much richer. We observed three different 

types of dynamics. In addition to coarsening, in the evolutionary relevant limit of 

C^ large times, oscillating traveling waves and heteroclinic orbits play a dominant role. 

J|T The weight of the different processes depends on the degree of mixing and the system 

^\ size. By analytical arguments and extensive numerical simulations we provide the full 

t" — characteristics of scenarios leading to extinction in one of the most surprising models 

'^ of ecology. 
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1. Introduction 

Stochastic many-particle systems provide a testing ground for non-equilibrium 
dynamics. In nature, systems frequently evolve away from equilibrium and then relax to 
an equilibrium steady state. Understanding the relaxation process is a central topic in 
non-equilibrium physics. Near-equilibrium fluctuations are governed by the same laws 
that hold in steady state and the transient is typically an exponential decay. Many 
systems, however, comprise absorbing states, which can be reached but never be left by 
the dynamics. In this case no fluctuations are present in the steady states. Such systems 
arise in a broad variety of problems, e.g. physics, chemistry or epidemics [T]. Much 
effort has been spent on the investigation of simple, diffusion-limited chemical reactions 
where the decay to equilibrium can obey power laws [2]. 

Understanding transitions into absorbing states is not only fundamental for non- 
equilibrium physics, but is also highly relevant for ecology. Here, absorbing states 
correspond to the extinction of species. Another characteristic feature of ecological 
systems are cyclic interactions. As a classic example, the work of Lotka and Volterra 
describes the dynamics of fish populations in the adriatic as persistent oscillations 
due to predator-prey interactions. Other examples include coral reef invertebrates [3j , 
rodents in the high arctic tundra in Greenland [4J, cyclic competition between different 
mating strategies of lizards [5] and chemical warfare of Escherichia coli bacteria under 
laboratory conditions [6]. 

Recent work has investigated cyclic competition in one-dimensional systems with no 
or only weak diffusion of the reacting agents [H El [9l [10]. Coarse-graining of temporally 
growing and annihilating domains has been identified as the mechanism that eventually 
leads to species extinction. However, individual's mobility may be significant and alter 
this picture qualitatively. 

In this article, we investigate the spatio-temporal dynamics of extinction in a 
paradigmatic model of three species in cyclic competition. Individuals are positioned 
on a one- dimensional lattice and are equipped with fast mobility that leads to effective 
diffusion. The system possesses absorbing states in the form of extinction of two of the 
three species and, because of fluctuations, the dynamics eventually comes to rest there. 
However, the time-scales until extinction occurs provide information on the stability 
of species diversity [llj. We identify three distinct types of dynamics that lead to 
extinction. These types of dynamics arise from the possible influences that intrinsic 
fluctuations can have on the coarsening process and on the traveling waves that the 
cyclic dynamics induces. The different dynamics lead to characteristic dependences 
of the extinction-time probability on the elapsed time t and the system size N. We 
provide semi-phenomenological arguments that quantify the functional form and the 
scaling behaviour of the extinction-time probability. These arguments yield information 
on the emergence and characteristics of the different types of dynamics. 
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2. The model 

Consider a stochastic, spatial variant of the May-Leonard model which serves as a 
prototype for cyclic, rock-paper-scissors-like species interactions. Three species A, B, C 
compete with each other in a cyclic manner, at rate a, and reproduce at rate fj, upon 
availability of empty space 0: 

AB A A0, BC A 50, CA A C0, 

Ads A AA, BdS A BB, CdS A CC. (1) 

For increasingly large populations intrinsic fluctuations eventually become 
negligible. If in addition spatial structure is absent, i.e., if every individual can interact 
with every other in the population at equal probability, the population dynamics is aptly 
described by deterministic rate equations for the densities s = (a, b, c) of the species A, B 
and C: 

dtSi = Si[fi{l- p) -aSi+2], for z e {1,2,3}. (2) 

Hereby the indices are understood as modulo 3 and p = a + b + c denotes the 
total density. May and Leonard showed that these equations possess 4 absorbing 
fixed points, corresponding to the survival of one of the species and to an empty 
system |12) . Furthermore a reactive fixed point s* = ^f^g (1, 1,1) exists that represents 
coexistence of all three species. Linear stability analysis shows that s* is unstable. The 
absorbing steady states that correspond to extinction, (1,0,0), (0,1,0) and (0,0,1), 
are heteroclinic points. The Lyapunov function C = abc/p^ demonstrates that the 
trajectories of the deterministic equations ([2]), when initially close to the reactive fixed 
point, spiral outward on an invariant manifold. On this manifold the trajectories then 
approach the boundary of the phase space and form heteroclinic cycles, converging to 
the boundary and the absorbing states without ever reaching them. 

However, intrinsic noise from finite-system sizes [T3| [H] and spatial correlations 
alter the above behaviour [T5| [TBI HTJ [TS]. While fluctuations ultimately drive the 
system into one of the absorbing fixed points [19], the formation of spatial patterns 
can substantially delay extinction and promote species coexistence [TBI [201 EI]- The 
resulting spatio-temporal dynamics of extinction is nontrivial and highly interesting. 

3. Numerical results 

We consider a one-dimensional lattice of L sites with periodic boundary conditions. Each 
lattice site hosts a fixed number M of individuals A, B, C and empty spaces 0, such that 
the concentrations in the rate equations ^ are given by the number of particles of a 
specified type divided by M. M may hence be viewed as the carrying capacity of a lattice 
site. The reactions ([l]) occur between individuals on the same lattice site. Individuals 
may change place with another individual or an empty space on a neighbouring lattice 
site at rate e. In order to keep the length scale, i.e. characteristic length scale for 
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Figure 1. (colour online) Spatio-temporal dynamics of three exemplary runs that 
correspond to the three classes of dynamics {M = 100, L ~ 600): rapid annihilation 
(left), heteroclinic orbits (center), and propagating waves (right). In the latter case an 
initially stable wave formation changes periodically at later times, leading to oscillating 
total densities. Colour encodes the concentrations of the three species (A, red; B, green; 
C, blue). 



diffusion, fixed when changing the lattice spacing L^^ we have to rescale e appropriately. 
In the continuum limit, where ^ holds, the exchange processes therefore lead to an 
effective diffusion of individuals at a diffusion constant D = eL^^ and thus to coupling 
between the lattice sites. In our simulations we implemented a continuous-time Markov 
process with sequential updating. At each simulation step an individual is randomly 
chosen. It then either reacts with a randomly chosen individual of the same site or 
changes place with a randomly chosen individual of the two neighbouring stacks, at 
probabilities corresponding to the rates a, fi and e. 

The population model introduced above possesses a net system size of A^ = ML 
which plays the role of an overall carrying capacity. For large enough M and L the 
intrinsic fluctuations have a strength proportional to the inverse square- root of N |22) . 
Different equivalent ways therefore exist for performing the thermodynamic limit, e.g., 
increasing the number M of individuals per lattice site, while keeping the lattice size L 
fixed or increasing the lattice size L, keeping M fixed. Because for large L a huge amount 
of exchange processes takes place between the reactions, requiring long computation 
time, we performed the thermodynamic limit in M — )■ oo and kept L = 100 fixed. The 
insensitivity of the results to the choice of the limit is supported by recent studies [lOj. 
L was chosen sufficiently large, such that L~^ was much smaller than the correlation 
length. 

We here consider the case of equal reproduction and selection rates fi = a = 1. 
Similar behaviour can be expected ioT fi ^ a, when D is rescaled appropriately |20| . and 
species dependent interaction rates fiU\ 123] . We chose a random initial configuration 
in which the density of the species approximately equals the ones of the internal fixed 
point s* of the rate equations ([2]). 

Our simulations reveal three distinct classes of dynamics (Figure II]). First, at 
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short time scales stochastic effects lead to the emergence of domains due to coarsening. 
In this scenario, after a short coarsening process, domains emerge, whose order does 
not correspond to the rules of cyclic dominance, leading to oppositely moving fronts 
and hence immediate annihilation. Extinction occurs rapidly in this scenario. The 
coarsening dynamics to extinction has been extensively studied [3 |B1 IS]. However, our 
simulations reveal a much richer dynamics to extinction. Two more processes dominate 
the dynamics for large times. Second, we observe situations where the population is 
almost entirely taken over by a species in cost of a second species, which dies out. 
A few individuals of the other surviving species are present in the system and, being 
the dominant one, slowly fixate. This scenario is intimately related to the heteroclinic 
orbits of the rate equations (pi). The global dynamics moves along the boundary of the 
invariant manifold of ([2]). Spatial patterns are of minor importance. Third, the system 
can enter a state of propagating waves of cyclically aligned, uniform domains. These 
states are only metastable: fluctuating front positions result in domain-annihilation and 
eventual extinction. For small D this effect was also denoted in |10] and corresponds to 
the spiral waves found in the two-dimensional model [20]. Rare events at the leading 
edge of the fronts cause the tunneling of domains and oscillating overall species densities. 
Figure |2(a)| provides a concrete picture of the different dynamical processes. 
The system's state probability is projected onto the invariant manifold of the rate 
equations ^. The colour signifles the logarithmic probability to flnd a given net density 
of species on the manifold before reaching an absorbing flxed point. We recognize that 
the system spends considerable time in the vicinity of the boundary, especially near 
the corners of the simplex, corresponding to the heteroclinic orbits occurring in the 
second scenario. The stochastic limit cycle around the unstable flxed point reflects the 
oscillating traveling waves from the third scenario (see also, e.g., |25]). In the center 
single trajectories of non-oscillating waves are visible. 



The influence of mobility is visualised in Figure 2(b) The Lyapunov function 
C = abc/p^ characterizes the system's behavior: it is zero at the boundaries and increases 
monotonically to the unstable flxed point s*. C therefore provides a measure for the 
distance of the system's state from the boundaries. The logarithmic probability to flnd 



the system at a certain value of £, depending on the diffusivity D, is given in Figure 2(b) 
A drastic change in the system's behavior occurs at a critical mobility Dc ^ 8 ■ 10~^. 
Above Df. we observe only heteroclinic orbits, characterized by a high probability to flnd 
the system at small values of C. For very small D the system exhibits traveling waves, 
performing random walks in concentration space. Below Dc both types of dynamics are 
present. The high probabilities for small values of C indicate heteroclinic orbits, while 
the ridge at larger values is caused by oscillating traveling waves. 

Quantiflcation of the three different dynamical scenarios is feasible through the 
extinction-time probability, P{t), meaning the probability density that two species go 
extinct at a certain time t. Mathematically it gives the probability distribution function 
of flrst-passage times into one of the absorbing flxed points. In our simulations we 
varied M from 1 to 2800 and set Z^ = 3 ■ 10~^, i.e. in the region, where all types 
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(a) 



(b) 



Figure 2. (colour online) | (a) [ Probability of net densities a, b, c for M = 600, L = 100, 
projected onto the invariant manifold of the rate equations (pi. Colour encodes the 
logarithm of the probability to find the system in the specific state, whereby red 
denotes the highest, yellow an intermediate, and blue a low probability. Note that 
the absorbing points themselves have not been included in the statistics. The graph 
allows to identify the reactive fixed point, an attractor for metastable oscillating waves, 
and the heteroclinic orbits. |(b)| The Lyapunov function C provides a measure for the 
distance of the system's state to the boundary of the simplex. The plot shows the 
logarithmic probability of net densities for different values of the diffusion constant 
D. Above a critical value of D we find heteroclinic orbits. For very small D we 
find traveling waves. Below Dc there is a region, where both heteroclinic orbits and 
traveling waves are present {M = 300, L = 100). 



of dynamics arise simultaneously. Figure 3(a) shows the extinction-time probability 



distribution for various system sizes. The sharp peak at small times results from the 
annihilation of oppositely propagating waves as a result of the coarsening process, i.e. 
the first scenario. The functional form of the extinction-time probability distribution for 
intermediate and large times is determined by the heteroclinic orbits and propagating 
waves, the second and third types of dynamics. We find an exponentially decaying tail 
that is, for large A^, preceded by a —3/2 intermediate asymptotic power-law interval. 
The length of this intermediate interval scales linearly with A^. The plateau or second 
maximum originates in the short term dynamics of the latter two scenarios. 

4. Semi-phenomenological arguments 



The characteristics of the critical behavior shown in Figure 2(b) can be understood 



through the spatial variant of the rate equations, ([2]). Following Ref. |26] the system's 
dynamics on the invariant manifold can, through a nonlinear transformation to variables 
za and zb, be recast in terms of the complex Ginzburg-Landau equation 



dtz = DV"^ + (ci - iuoQ)z - C2(l + ic^] 



\z?z. 



(3) 
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Figure 3. (colour online) (a) Double-logarithmic plot of the extinction-time 
distribution P(t) for several system sizes. A sharp peak at small times is followed 
by a second maximum or plateau and by an intermediate t~^'^ power law. The length 
of th e power law region scales with N. The tail of the distribution decays exponentially. 



(b) Semi-logarithmic plots of the survival probability S{t) = 1 — Jp P{t')dt' for different 
N . S{t) exhibits the same long-time exponential decay as P{t). With the rescaling 
t/N for large (top right, N = 30000 to iV = 280000) and t\n{N)-^ for small systems 
(bottom right, N — 5000 to A^ = 20000) the exponential tails collapse onto universal 
curves, in agreement with our analytical predictions. 



with ci 



fia 



2(3M+ff) ' 



C2 = 



o-(3m+o-)(48^+11o-) 
56tJ.{3fi+2a) 



and C3 



ASfi+Ua 



|27| . The theory of front 



propagation into unstable states predicts that (pi) always admits traveling waves as stable 
solutions |28]. Following a classic treatment of the problem of front-speed selection 



we obtain their wavelength as A 



2Trc:i^/D 



At the critical diffusivity Dc the 



wavelength A exceeds the system size such that the fronts become unstable. From the 
condition A = 1 and accounting for the rescaling factor mentioned in [27] we obtain 
Dc ~ 7.6 ■ 10"'', which is in very good agreement with our numerical results. 

The behaviour of the extinction-time probability distribution can be understood 
through semi-phenomenological models. In the following we show how such models 
yield the shape of the extinction-time probability distribution and its dependence on 
A^. In particular, we give an explanation for the scaling behaviour of the power-law 
interval and the long-time exponential decay. We show that, depending on the system 
size, either heteroclinic orbits or traveling waves dominate the long-time dynamics. 

The dynamics of the heteroclinic orbits can be quantified as follows. Consider a 
small density Oq G 0{N~^) of individuals of species A in a large pool of species B 
[b G 0(1)]- Due to reproduction of B empty space is as sparse as A individuals are: 
1 — a — b E 0{N~^). Simulations inform us that spatial patterns are not relevant in this 
scenario. We therefore consider a well-mixed system of size A^. Three reactions lead to 
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the take-over of the population through the dominating species A: 

AB ^ A0 at rate Naah G C(l), 

50 ^ BB at rate N ^ih{l - a - b) e 0(1), 

A0 ^ AA at rate Nna{l - a - b) E C(iV"^). 

Here the rates are meant as transitions per unit time. The fast processes are in 
equihbrium and can be adiabatically ehminated for large N, yielding a = -{1 — a — b). 
Species A occurs at the same density as empty sites times fi/a. A series of reproduction 
events remains, each with an exponentially distributed waiting time. In the language 
of stochastic processes this is a pure birth process, studied arising in preferential 
attachment problems. The extinction-time probability distribution Ph{t\ao) is thus 
given by a convolution over exponential functions. By applying the Laplace transform 
one can show that it can be expressed in closed form as 

N{l-ao) N{l-ao) 

P,{t;ao)= Y. ^^^~''' n v^' (4) 

i=0 i,j=0,i¥'j 

with rates A,, = Na (ao + i/N) |29j. One finally has to marginalize over the probability 
p(ao) of starting with an initial density oq of species A to obtain the extinction-time 
probability distribution Ph{t) as results from the heteroclinic-orbit dynamics: 

P,{t) = Y.Ph (t; «o = ^) P (ao = ^) . (5) 

For any reasonable p(ao) the asymptotic behaviour is dominated by the term for the 
lowest initial concentration ao = 1/A^ and the lowest reproduction rate Aq: Ph{t) ~ 
exp (—at/N) , for t — > oo. We thus find a A^~^-dependence of the exponential tail 
on the system size. For intermediate times we have to take the full convolution and 
marginalization sums of eq. (IS]) into account. Numerical evaluation indeed yields a 
—3/2 power-law interval for uniformly distributed ag. The time of crossover between 
the power-law and the exponential decay scales with A^. We therefore find that the 
second type of dynamics, heteroclinic orbits, lead to the intermediate power-law regime 
in the extinction-time probability distribution. The exponential tail of this distribution 
can either result from heteroclinic orbits or from propagating waves as shown below. 

The third type of spatio-temporal dynamics, propagating waves, is metastable. 
They disappear only through the rare annihilation of neighbouring fronts. By symmetry, 
the waves move with the same average velocity. For small D the domain interfaces are 
sharp and therefore interact only on distances that are much smaller than the average 
domain size. The extinction dynamics in this scenario can therefore be described 
within an interface picture, where, in a comoving frame, the wave fronts behave 
as random walkers on a one dimensional lattice with diffusion coefficient Df. For 
larger D long-range interactions between the interfaces become important, leading to 
a tunneling of domains. However, within the interface picture this merely corresponds 
to a relabeling of interfaces and therefore does not influence the extinction dynamics. 
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Numerical simulations validate the assumption of normal diffusion. For a single series 
of subsequent A,B,C domains the survival probability Su,{t), meaning the probability 
that all three domains still coexist at time t, follows as the survival probability of a single 
random walker between absorbing boundaries at distance /. The probability distribution 
Cw{t, x; Xq, I) for the random walker to be at position x and time t when starting at Xq 
obeys a diffusion equation subject to absorbing boundary conditions. The solution is 
well known: 

°° 
cUt, X- xo, l) = Y,^n sin (^) e-^^) ^^* , (6) 

n=l 

with the coefficients An = |;sin(^^^^) being determined by the initial condition 
Cw{x, t] Xq, I) = 6{x — Xq), scc e.g. j30]. Averaging over space, the initial positions xq and 
identically distributed interval lengths / yields the survival probability S^it) from the 
traveling-wave dynamics: 

„ oo ^ „i 

In the asymptotic limit the expression evaluates to S'^(t) ~ e~^^f'" *, for t — > oo. 
The extinction-time probability distribution follows as Pu]{t) = —dSw{t)/dt. What is 
the diffusion constant Df of the domain front? Brunet et al. proposed [HI] that the 
diffusion constant for a broad class of stochastic propagating waves depends on A^ as 
Df ~ In(A^)^'^. Therewith the exponential decay of the extinction-time probability 
distribution's tail, as resulting from the traveling-wave dynamics, is proportional to 



ln(A^)~^. This result is validated by our numerical simulations, see Figure 3(b) bottom 
right. Heteroclinic orbits and traveling waves both contribute to the asymptotic limit of 
the net extinction-time probability distribution P{t): P{t) ~ P/i(t)-|-P^(t), for t — )■ oo. 
Both contributions yield exponential decays at large times, but with different scalings 
in A^. For small systems the ln(A^)~^ term in P«,(t), resulting from the traveling- wave 
dynamics, dominates. In contrast, the 1/A^-decay in Ph{t) resulting from heteroclinic 
orbits yields the major contribution when A^ is large. In agreement with these analytical 
results we indeed find numerically that the exponential decay scales as ln(A^)~^ for small 



A^ and as 1/A^ for large A^ (Figure 3(b) ). Numerically we identified the crossover between 



both regimes to occur aX N '^ 20, 000. 

5. Conclusion 

We investigated the spatio-temporal extinction dynamics in a three species stochastic 
population model with cyclic interactions. While previous work has mainly focused 
on the coarse graining process that drives the system to extinction we identified two 
more types of dynamics that are rare but due to their lifetime most important from 
an evolutionary perspective. The three classes of dynamics, namely rapid annihilation 
of domains, heteroclinic orbits, and traveling waves are correlated with features of the 
phase portrait and leave their fingerprints in the extinction-time probability distribution. 
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The weight of these processes depends on the degree of mixing as well as on the system 
size. Based on the different dynamical scenarios we provided semi-phenomenological 
calculations that yield the functional form of this probability distribution and its 
dependence on the system size. We believe that our results are of general relevance as we 
expect a similar phenomenology in other systems described by the complex Ginzburg- 
Landau equation. 
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